{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predicting Porosity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example explains how to estimate porosity of a network. We also discuss some challenges in estimating the porosity of the network and how to reduce the estimation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import openpnm as op\n",
    "import porespy as ps\n",
    "import matplotlib.pyplot as plt\n",
    "ps.settings.tqdm[\"disable\"] = True\n",
    "op.visualization.set_mpl_style()\n",
    "np.random.seed(10)\n",
    "%matplotlib inline\n",
    "np.set_printoptions(precision=5)\n",
    "op.Workspace().settings.loglevel = 40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a random cubic network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pn = op.network.Cubic(shape=[15, 15, 15], spacing=1e-6)\n",
    "pn.add_model_collection(op.models.collections.geometry.cubes_and_cuboids)\n",
    "pn.regenerate_models()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate void and bulk volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of Porosity is: 0.21\n"
     ]
    }
   ],
   "source": [
    "Vol_void = np.sum(pn['pore.volume'])+np.sum(pn['throat.volume'])\n",
    "inlet = pn.pores('left')\n",
    "outlet = pn.pores('right')\n",
    "A = op.topotools.get_domain_area(pn, inlets=inlet, outlets=outlet)\n",
    "L = op.topotools.get_domain_length(pn, inlets=inlet, outlets=outlet)\n",
    "Vol_bulk = A * L\n",
    "Poro = Vol_void / Vol_bulk\n",
    "print(f'The value of Porosity is: {Poro:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussions and Issues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Domain volume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the issues in estimation of porosity of the network is to estimate the domain volume correctly. In a cubic network for example, finding the length at x direction using Nx * spacing is erroneous. The reason is the domain length in x direction includes additional lengths from half of pore diameter for pores that locate on the left and right side. This is shown in figure below:\n",
    "\n",
    "left) the green plane is located at the pore centers of left boundary pores and some pore volume on the left is ignored. Applying a similar plane on other sides of the network to find the length in each direction, the resulting bulk volume could be underestimated. \n",
    "\n",
    "right) the green plane is located at the far left side. Applying a similar plane on other sides of the network to find the length in each direction, the resulting bulk volume could be overestimated. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image](https://user-images.githubusercontent.com/43128873/187976526-d36c60be-8514-4326-8ee5-00a24ea90d56.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overlapping pores and throats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another issue is to ensure that the pore-scale models for volume are consistent. For example there is no overlapping pores, because in this case the void volume calculation will overestimate the real void space. Depending on the pore scale model used for the shape of pores and throats, they may need special methods to calculate the void volume to account for the overlap between throats and pores. Depending on the method that was used for assigning throat lengths, this overlap volume may be included in volume calculations. Existing methods to correct the throat volumes are ``lens`` and ``pendular_ring`` methods available in geometry models for throat_volume. For example assuming a spheres and cylinders geometry, `lens` model in geometry collection tackles this problem. The `spheres_and_cylinders` geometry collection includes `throat.total_volume` and `throat.volume`. `throat.volume` is the volume of throat with corrections using ``lens`` volume:\n",
    "\n",
    "Throat volume (`throat.volume`) = volume of cylinder (`throat.total volume`) - the overlap of cylinder with spherical pores at its two ends  (`difference`)\n",
    "\n",
    "Let's create a spheres and cylinders geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)\n",
    "pn.regenerate_models()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial Porosity: 0.12852\n",
      "Corrected Porosity: 0.12590\n"
     ]
    }
   ],
   "source": [
    "Vol_void_initial = np.sum(pn['pore.volume'])+np.sum(pn['throat.total_volume'])\n",
    "Vol_void_corrected = np.sum(pn['pore.volume'])+np.sum(pn['throat.volume'])\n",
    "Poro_initial = Vol_void_initial / Vol_bulk\n",
    "Poro_corrected = Vol_void_corrected / Vol_bulk\n",
    "print(f'Initial Porosity: {Poro_initial:.5f}')\n",
    "print(f'Corrected Porosity: {Poro_corrected:.5f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although in this example the lens volume was negligible, depending on the size of the pores and throats, this value can be too high to be neglected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. Note:: <Pendular ring correction>\n",
    "\n",
    "    `Pendular_ring` method calculates the volume of the pendular rings residing between the end of a cylindrical throat and spherical pores that are in contact but not overlapping. This volume should be added to the throat volume if the throat length was found as the center-to-center distance less the pore radii."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracted networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Void volume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In extracted networks, different geometrical shapes can be assigned to the pores and throats and their volume can be calculated using existing geometry models. However, the original segmented pore regions are not regular shapes. Therefore, the total pore and throat volumes in the network don't add up to the known void volume, since the pores and throats don't fill space perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity from image: 63.3%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAOgCAYAAAAjxRuxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAB7CAAAewgFu0HU+AABNhklEQVR4nO3dfZCU1Z0v8N8AAwgoOrwJEYjEgOwSI0NUlBIQJCRGLSVZMJqrrC+oa625WtcETTQS3bir2VRysxVj4i0grkGMUZPVJERBRVFMFFnRBWJQg8hrRJGXGRih7x8UT2Zghmlmus+8fT5VVJ2n+zxnTjNnuuc7p/v5leRyuVwAAABAkbVr6gkAAADQNgigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASLT6A7tq1K+67774466yzYuDAgdG5c+fo27dvnHbaafHd7343/vrXvzb1FAEAAIiIklwul2vqSTTUihUr4sILL4xXXnmlzj69e/eOmTNnxllnnZVwZgAAAOyvxQbQNWvWxCmnnBJr166NiIiSkpIYPXp0HHfccbFx48Z48skno6KiIiIiSktL47e//W2MHz++KacMAADQprXYADpmzJhYuHBhREQMHDgwfv3rX8cJJ5yQ3f/Xv/41Lrjggpg/f35ERJSVlcWqVaviyCOPbIrpAgAAtHkt8jOgv/nNb7Lw2bFjx/iv//qvGuEzIqJnz57xq1/9KgYNGhQREZs3b44777wz+VwBAADYq0XugH7hC1+I3/zmNxERccUVV8RPfvKTOvvef//98ZWvfCUi9u6CbtiwITp06HDQ8Y8//vh49913a9zWpUuXLMwCAAC0Zm+++Wbs2LGjxm0f+9jHYsWKFY0at8UF0G3btkXPnj1j586dERHx/PPPx6mnnlpn/507d0avXr1i69atERExf/78GDdu3EG/xuGHHx7btm0r3KQBAABauG7dumW5qqFa3Ftwn3/++Sx8du3aNU466aSD9u/UqVOMHDkyO16wYEFR5wcAAEDtWlwAXb58edb+1Kc+Ve/baSMiysvLaz0fAACAdFpcAF25cmXWHjhwYF7nDBgwIGs39j3LAAAANEz924fNzHvvvZe1+/Tpk9c5Rx99dNbevHlzvf27dOlS0M+Ado+ygo3VUENHbG/qKVAkr63YFdu21/wod7euJTHs+I5Jvv7yl7sm+TqtQVv7OWzqtQl1sTZpzqxPCm3xyzsLNlaXLl0aPUaLC6DVg+Fhhx2W1znV++UTLAcNGhQbN2489MnV4aSSg1/0KIV5jy1t6ilQJKPOfueAJ5Zhx3eMRY/1T/L1J/Y7McnXaQ3a2s9hU69NqIu1SXNmfVJo7fv+uWBjFaIqSIt7C25lZWXW7tgxv78EderUKWtXVFQUfE4AAADUr8UF0M6dO2ftXbt25XXOvqvmRuS/awoAAEBhtbi34Hbr1i1r57ubWb1f9fMbo320j27RvSBjAQAAFMPIEZ3q71SL2j6PXAgtLoD26NEja2/YsCGvc9avX5+1y8oKc0GgbtG9WXy2M1/5fE5v3tqlRZ8HLY/PeBaOn8OWJ8X69z0HoJga+vnh2j6PXAgt7i24Q4YMydp/+ctf8jpn9erVWfv4448v+JwAAACoX4sLoEOHDs3ay5Yti48++qjec5YsWVLr+QAAAKTT4gLoaaedll3Vdvv27fHSSy8dtP/OnTtj8eLF2fG4cS3nbbMAAACtSYsLoN26dYvx48dnx7NmzTpo/4cffji2bt0aERFHHXVUjB49upjTAwAAoA4tLoBGRPzTP/1T1p45c2a8/vrrtfbbsWNH3HLLLdnxlVdeGR06tLjrLgEAALQKLTKAfuELX4jTTz89IvbWAj377LNj2bJlNfq89957cd5558Wf//zniNh79duvf/3ryecKAADAXi12O/DnP/95nHzyybFu3bp4++2348QTT4wxY8bEoEGDYtOmTfHkk0/Gjh07IiKiQ4cO8eCDD8aRRx7ZtJMGAABow1psAD3mmGNiwYIF8eUvfzmWLl0ae/bsiaeeeiqeeuqpGv169eoVM2fOrPG5UQAAANJrsQE0Ym9NzxdffDEeeOCBmDNnTrz++uuxYcOGOPLII2PQoEFx/vnnx6WXXho9e/Zs6qlCszSx34lNPQX2U9/3ZN7apUnm0RK0lvWb4nFYNwA0Fy06gEZEdOzYMS6++OK4+OKLi/p1ukdZnFSihAvNz6LH+jf1FKBW1ibNlbVJc2Z90tq1yIsQAQAA0PIIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAk0eLLsFA46g+2Pq2lTiJ/05Z+Tq3fwsnn/7I1rR0Ami87oAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEuqAQgulRiK1aSn1Hq3f5qct1ZgFoOnYAQUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASKJDU08AgNZnYr8Tm3oKFFh939N5a5cmmQcALZsdUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCXVAATgkanxSm3zWhVqhANgBBQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIokNTTwCoXT5F3aEh6ltb89YuTTIPAKDtsQMKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJNGhqScAQPMysd+JTT0FWqn61ta8tUuTzAOApmMHFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgiQ5NPQGaDwXAm5f6vh/1FXSHulhbAEBTsQMKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJNGhqSdAGvPWLm3qKQAtRH3PFxP7nZhkHrQ+XosAsAMKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJBEh6aeAPVTuJva5LMuJvY7sejzoOWxLgCApmIHFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQh3QPA0dsT3mPba0qacBh6S+WqHqQQIAkJIdUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJDo09QSApjNv7dJ6+0zsd2LR5wG0fPk8nwCAHVAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAl1QKENU+MTKJR8nk/UCgXADigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJJAmgb7/9dvz0pz+Nr3zlK/HpT386jjrqqCgtLY2ysrI44YQT4sorr4xnnnmmQWPPnz8/Lr744hg8eHB07do1G/OGG26IFStWFPiRAAAA0FAdijn4K6+8EldddVX84Q9/qPX+999/P95///1YtmxZ/OQnP4mxY8fG7NmzY8CAAfWO/eGHH8a0adNi7ty5NW7fsWNHNuYPfvCDmDFjRtx4440FeTwAAAA0XFED6MqVKw8In4MHD45hw4ZFz54944MPPojnn38+1qxZExERTz/9dJx66qnx7LPPxqBBg+oct6qqKiZNmhTz58/Pbhs2bFiMGDEiKioqYuHChbF+/fqoqqqKm266KaqqquKWW24pzoMEAAAgL0UNoPscd9xxcfnll8dXvvKV+NjHPlbjvj179sTMmTPj2muvjR07dsTatWvjoosuiueffz5KSkpqHe+2227Lwmfnzp1j5syZccEFF2T379q1K775zW/GXXfdFRER3/rWt2LMmDExZsyYIj1CAAAA6lPUz4D27ds3Zs6cGStWrIivf/3rB4TPiIh27drFZZddFv/5n/+Z3bZ48eL4/e9/X+uYGzdujO9973vZ8fe///0a4TMiomPHjnHnnXfGlClTstu8DRcAAKBpFTWAjhkzJqZOnRrt27evt+/5558fJ598cnb8+OOP19pv9uzZsX379ojY+3beadOm1TnmnXfeGe3a7X2IL7zwQrzyyiuHMn0AAAAKqFmVYRk1alTWfvvtt2vt8+ijj2btqVOn1vk23YiIAQMGxPjx47PjRx55pNFzBAAAoGGaVQCtHiZ37959wP2VlZWxePHi7Hjs2LH1jlm9z4IFCxo1PwAAABquWQXQZcuWZe3+/fsfcP/KlStjz549EbE3rA4fPrzeMcvLy7P28uXLCzBLAAAAGiLJVXDz8c4779TYoTzzzDMP6LNy5cqs3bt37+jcuXO941avKbp58+bYtGlT9OrV65Dn99qKXTHq7HcO+byIiEWPHRimAQAAiq2hGea1FbsKPJO9mk0Ave6667K33Q4YMCDOOeecA/q89957WbtPnz55jXv00UfXON68eXODAui27blY/PLOQz4PAACgqTS3DNMs3oI7e/bs+OUvf5kd33HHHdGpU6cD+m3bti1rH3bYYXmNvX+/6mMAAACQTpMH0Jdeeimuuuqq7HjKlClx4YUX1tq3srIya3fs2DGv8fcPshUVFQ2YJQAAAI3VpAH0rbfeinPOOScLlp/61KfinnvuqbN/9c987tqV33uSd+6sueWc784pAAAAhdVknwFdt25dTJgwIdavXx8REYMGDYp58+ZF9+7d6zynW7duWTvfncz9+1Uf41B061oSw47Pb9cVAACgORg54sCPNubjtRW7Ytv2XIFn00QB9L333osJEybEqlWrIiKib9++8eSTT0bfvn0Pel6PHj2y9oYNG/L6WvsC7j5lZWWHONu9hh3f0dVsAQCAFqWhGWbU2e8U5QJGyd+C++GHH8bnPve5eP311yNib6h84okn4thjj6333CFDhmTtjRs31vhMaF1Wr16dtcvKyhp0BVwAAAAaL2kA3b59e5x11lnx0ksvRUTEEUccEfPmzYu///u/z+v8IUOGRLt2e6ecy+Vi6dKl9Z6zZMmSrD106NBDnzQAAAAFkSyAVlZWxrnnnhuLFi2KiIguXbrEb37zmxgxYkTeY3Tu3DlGjhyZHT/99NP1nvPMM89k7XHjxuU/YQAAAAoqSQCtqqqKL37xi7FgwYKI2Fsa5Ve/+lWMGjXqkMc677zzsvasWbMO2nfNmjUxf/78Ws8FAAAgraIH0N27d8eFF14Yv/nNbyIiokOHDvHggw/GmWee2aDxLrnkkujatWtERKxcuTLuvffeOvt+7Wtfi927d0dExKmnnhrl5eUN+poAAAA0XlEDaC6Xi8svvzweeuihvV+sXbu477774txzz23wmL17947rr78+O7722mvjwQcfrNFn165dMX369JgzZ0522x133NHgrwkAAEDjFbUMy913313jbbKf+MQn4rnnnovnnnuu3nN79OgRM2bMqPW+m2++ORYtWhQLFiyIioqKmDJlStx+++1RXl4elZWVsXDhwli3bl3Wf8aMGTFmzJhGPx4AAAAarqgBdOPGjTWO33jjjXjjjTfyOnfgwIF1BtDS0tJ4+OGHY9q0adnu57Jly2LZsmUH9Lv11lvjpptuasDsAQAAKKSiBtBi6t69e8ydOzeuuOKKmD17drzwwguxbt26KC0tjf79+8fEiRPjsssuU3oFAACgmShqAL311lvj1ltvLeaXiDPPPLPBFzQCAAAgnWR1QAEAAGjbBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkmixdUABaLnmrV1ab5+J/U4s+jwonHy+pwBgBxQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgiQ5NPQGg6cxbu7TePhP7nVj0edD65LO2GjuGtZlWIb6nAGAHFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgiQ5NPQGgeauv+PzEficmmQfNR31rIpVCzKOtrN/m8j0DADugAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEAS6oAC0Gapj0lzlaJGrfUPNAU7oAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEuqAAo1SiDpyKerd8Tdq/0FxtZTntELM0/MJcKjsgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJNGhqScAUF8h85ZS1L05UBQeGs9zTv4K8X/leQvaFjugAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEAS6oDSoqSqzaYmWfPi+wHQdqV47fc6A+nYAQUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAklAHlKRS1fFsrMbOUz0xgOappbwOtSXN4XtSiDl47Yf82AEFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEiiQ1NPgJajORSKbikUtAY4dF5nAFo/O6AAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBLqgAIASajzSWtW3/puS/W9U/yst6X/z9bGDigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBIdmnoCNB8KhDcvCloDLYnXEGgdWsrPciHm6XeppmEHFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQh1QAKBeLaU2IADNmx1QAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJdUBbALXXqE0+62Le2qVFnwfQOnitgcZpKa+5ftb/pr7/i5byPW1p7IACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACTRoakn0BYo+AsAAGAHFAAAgEQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQh3QRlLjEwAAID92QAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJNQBhRZq3tqlTT0FoIVQsxrYp77fH9rS84XfpZqGHVAAAACSaPIAet1110VJSUn27+Mf//ghnT9//vy4+OKLY/DgwdG1a9coKyuLE044IW644YZYsWJFcSYNAADAIWvSt+D+4Q9/iP/7f/9vg8798MMPY9q0aTF37twat+/YsSPef//9WLZsWfzgBz+IGTNmxI033liI6QIAANAITRZAq6qq4vLLL489e/Y06NxJkybF/Pnzs9uGDRsWI0aMiIqKili4cGGsX78+qqqq4qabboqqqqq45ZZbCjl9AAAADlGTvQX33/7t32LZsmUREXHhhRce0rm33XZbFj47d+4cc+bMiWXLlsWsWbNi7ty58Ze//CVuuOGGrP+3vvWteOaZZwo3eQAAAA5ZkwTQFStWxO233x4RERdddFFMmDAh73M3btwY3/ve97Lj73//+3HBBRfU6NOxY8e48847Y8qUKdlt3oYLAADQtJIH0FwuF5dffnns3LkzjjrqqBphMh+zZ8+O7du3R0TE4MGDY9q0aXX2vfPOO6Ndu70P8YUXXohXXnml4RMHAACgUZIH0LvvvjsWLVoUERF33XVX9O7d+5DOf/TRR7P21KlTo6SkpM6+AwYMiPHjx2fHjzzyyKFNFgAAgIJJehGiNWvWxPTp0yMi4vTTT49LL730kM6vrKyMxYsXZ8djx46t95yxY8fGE088ERERCxYsiG9/+9uH9DX3Wf5y1zZVmBcAgPzV93vivLVLk8wDmrukO6BXX311bN26NTp27Bj33HPPQXcva7Ny5crsqrklJSUxfPjwes8pLy/P2suXLz+0CQMAAFAwyXZAH3jggXjsscciIuLrX/96DB069JDHWLlyZdbu3bt3dO7cud5zBgwYkLU3b94cmzZtil69eh3y194WW+KPuQWHfF5ExEkl4xp0HgAAQGOMOvudBp332opdBZ7JXkkC6HvvvRdf/epXIyLik5/8ZHzjG99o8Dj79OnTJ69zjj766BrHmzdvblAA3R27Y0tsPuTzAAAAmsril3c29RRqSPIW3Ouuuy42btwYERH33HNPdOrUqUHjbNu2LWsfdthheZ2zf7/qYwAAAJBO0QPo73//+7jvvvsiIuKSSy6JM844o8FjVVZWZu2OHTvmdc7+YbeioqLBXx8AAICGK2oA3b59e1x55ZUREdGjR4/47ne/26jxqn/mc9eu/N6TvHNnzS3nfHdOAQAAKKyifgb0G9/4Rrz99tsREfHv//7v0bNnz0aN161bt6yd707m/v2qj3Eo2kf76BbdG3QuAABAUxg5omEff3xtxa7Ytj1X4NkUMYAuWbIkfvjDH0ZExBlnnBGXXHJJo8fs0aNH1t6wYUNe56xfv77GcVlZWYO+drfo7mq2JKVeGAC0HvnUk28Or/2FmEM+j7WxmsP/VUux6LH+DTpv1NnvFOUCRkULoK+++mpWs3P16tUxcuTIOvtu2rQpa69bt65G35tvvjm+8IUvRETEkCFDsts3btwYlZWV9ZZiWb16ddYuKytr0BVwAQAAaLwkZVhWrVoVq1atyqvvrl274sUXX8yOq4fTIUOGRLt27WLPnj2Ry+Vi6dKlBw22EXt3YvdpSO1RAAAACiNJGZZC6dy5c43A+fTTT9d7zjPPPJO1x43zFloAAICmUrQAOnXq1Mjlcnn9mzlzZnbewIEDa9w3derUGuOed955WXvWrFkHncOaNWti/vz5tZ4LAABAWi1qBzRiby3Rrl27RkTEypUr4957762z79e+9rXYvXt3RESceuqpUV5enmSOAAAAHKjFBdDevXvH9ddfnx1fe+218eCDD9bos2vXrpg+fXrMmTMnu+2OO+5INkcAAAAOlOQiRIV28803x6JFi2LBggVRUVERU6ZMidtvvz3Ky8ujsrIyFi5cGOvWrcv6z5gxI8aMGdOEMwYAAKBFBtDS0tJ4+OGHY9q0adnu57Jly2LZsmUH9Lv11lvjpptuaoppAgAAUE2LDKAREd27d4+5c+fGFVdcEbNnz44XXngh1q1bF6WlpdG/f/+YOHFiXHbZZUqvANDm5VOwPUXheKBtyOc5h7arWQTQqVOnHnC123ydeeaZceaZZxZ2QgAAABRci7sIEQAAAC2TAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEs2iDAu0NupfAQDAgeyAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJqAMKtVDHEwAACs8OKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEh2aegIAQNObt3bpQe+f2O/EJPOAtqyt/JzV93xD62YHFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQh1QAKBe6oQChVKI5wu1RFsuO6AAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBLqgOZp6IjtMe+xpQfcru5Zy6NuFEDh5fPc6jUTKJT6nk/8vtd82QEFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEiiQ1NPAABoG+orDF9fYXmAfOXzfFLfcxLFYQcUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCHdBGyqd+kLpmaanpBNAypXj+9poM0LTsgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASagDmkBj65q1pZplangCUEyFeJ1pS6/L0JrV97Ps99LisAMKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJBEh6aeAPVTBBcAmo/6XpfrK24P0JbZAQUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIokNTTwAAoDWZt3Zpo8eY2O/ERo8B0BzZAQUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAklAHFACgmamvlqg6oUBLZQcUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCHVAAgBamvjqh1KRuKjQfdkABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJLo0NQTAIpH4e38KeoOAFB8dkABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCTUAYVmSg3PtArx/62WKEDzVN/zs9fctsnrdtOwAwoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJKEOKBSBemJtU33fd/XGAIC2zg4oAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASHZp6AgBtxcR+J9bbZ97apUWfBxRDPuu7tfBz2vrU9z1tS+u7tfBz2nwl3wFdsmRJTJ8+PT7zmc9E3759o1OnTtGvX78oLy+PSy+9NO67775Yv359XmPNnz8/Lr744hg8eHB07do1ysrK4oQTTogbbrghVqxYUeRHAgAAwKFItgO6cePGuP766+P+++8/4L5169bFunXr4pVXXomZM2fGNddcE//xH/9R51gffvhhTJs2LebOnVvj9h07dsT7778fy5Ytix/84AcxY8aMuPHGGwv+WAAAADh0SQLo6tWrY+zYsfHWW29ltx177LFRXl4ePXr0iIqKinjjjTdi6dKlUVlZedCxqqqqYtKkSTF//vzstmHDhsWIESOioqIiFi5cGOvXr4+qqqq46aaboqqqKm655ZaiPTYAAADyU/QAumXLljjjjDOy8FleXh4//OEP47TTTjug77Zt2+Lxxx+PXC5X53i33XZbFj47d+4cM2fOjAsuuCC7f9euXfHNb34z7rrrroiI+Na3vhVjxoyJMWPGFPJhAQAAcIiKHkD/z//5P/Hmm29GRMTo0aPjt7/9bXTp0qXWvt26dYspU6bUOdbGjRvje9/7Xnb8/e9/v0b4jIjo2LFj3HnnnbF69ersLbo33nhjPP/88419KAAAADRCUS9CtHTp0rj33nsjIuLwww+P+++/v87wmY/Zs2fH9u3bIyJi8ODBMW3atDr73nnnndGu3d6H98ILL8Qrr7zS4K8LAABA4xU1gP74xz/O2pdeemkcc8wxjRrv0UcfzdpTp06NkpKSOvsOGDAgxo8fnx0/8sgjjfraAAAANE7R3oK7e/fumDNnTnZ80UUXNWq8ysrKWLx4cXY8duzYes8ZO3ZsPPHEExERsWDBgvj2t7/dqDnAPuqBUSz1rS11zWgqnvf+xs9p25PP99TPSFp+zlquou2Avvbaa/Hhhx9GRETXrl1j+PDhsXPnzrjnnntizJgx0bt37+jcuXMcc8wxcfbZZ8dPf/rT2LVrV53jrVy5Mvbs2RMRESUlJTF8+PB651BeXp61ly9f3shHBAAAQGMUbQf0j3/8Y9YeMmRIrFq1Kr70pS/Fa6+9VqPfu+++G++++248/vjjcccdd8RDDz1UIzjus3Llyqy9L7zWZ8CAAVl78+bNsWnTpujVq1dDHk68tmJXjDr7nQadu+ix/g06DwAAoDEammFeW1H35mBjFC2AvvPO3x5ou3bt4rOf/WysXr06IiKOP/74OOmkk6J9+/bx6quvxpIlSyIi4q233orRo0fHs88+e8AO53vvvZe1+/Tpk9ccjj766BrHmzdvbnAA3bY9F4tf3tmgcwEAAJpCc8swRQugH3zwQdZ+6aWXIiLisMMOi1mzZsXkyZNr9H3qqadi8uTJ8de//jW2b98eU6ZMiddffz1KS0uzPtu2bcvahx12WF5z2L9f9TEAAABIq2ifAd1XLqW62bNnHxA+IyLOOOOM+PWvf52VTXnjjTfi/vvvr9GnsrIya3fs2DGvOXTq1KnGcUVFRV7nAQAAUHhFC6D7f0bzpJNOin/4h3+os/+pp54akyZNyo4feOCBOsc72MWKqtu5s+Z2c747pwAAABRe0d6C261btxrH559/fr3nnH/++fHQQw9FRMTzzz9f53j57mTu32//OR2Kbl1LYtjx+e28AgAANAcjR3Sqv1MtXluxK7ZtzxV4NkUMoD169Khx/Hd/93f1nlO9z9atW2Pr1q1x+OGHHzDehg0b8prD+vXraxyXlZXldV5thh3f0dVsAQCAFqWhGWbU2e8U5QJGRQugxx9/fI3jfHYf9+9TPYAOGTIku33jxo1RWVlZbymWfVfdjdgbPht6BVwAaO0m9juxqafQauTzfzlv7dKiz4O0CvE9bSs/h9Z/21a0z4AOGzasxvHWrVvrPWf/Pt27d8/aQ4YMyS5SlMvlYunSpfWOt6+8S0TE0KFD6+0PAABA8RQtgB577LExaNCg7Ph//ud/6j2nep+ysrLo2rVrdty5c+cYOXJkdvz000/XO94zzzyTtceNG1dvfwAAAIqnaAE0ouaFhx555JF6+1fvM3r06APuP++887L2rFmzDjrWmjVrYv78+bWeCwAAQHpFDaBXX311lJaWRkTESy+9FL/4xS/q7PvCCy/UCKBTp049oM8ll1yS7YquXLky7r333jrH+9rXvha7d++OiL0lXsrLyxvyEAAAACiQogbQT3ziE/FP//RP2fHUqVNrDaFPPfVUnHvuubFnz56IiBg5cmSce+65B/Tr3bt3XH/99dnxtddeGw8++GCNPrt27Yrp06fHnDlzstvuuOOORj8WAAAAGqdoV8Hd59/+7d9iyZIl8eyzz8aOHTti8uTJMXTo0DjppJOiffv28eqrr8bLL7+c9e/bt288+OCDUVJSUut4N998cyxatCgWLFgQFRUVMWXKlLj99tujvLw8KisrY+HChbFu3bqs/4wZM2LMmDHFfpgAAADUo+gBtFOnTvFf//VfcfXVV2e7ksuXL4/ly5cf0PeUU06JX/ziF9G/f921akpLS+Phhx+OadOmZbufy5Yti2XLlh3Q79Zbb42bbrqpgI8GAACAhip6AI3YW07l5z//eVx11VXxs5/9LJ577rl49913Y/fu3dGnT58YOXJkTJ48Oc4777w6dz73H2/u3LlxxRVXxOzZs+OFF16IdevWRWlpafTv3z8mTpwYl112mdIrNFhbqcNFy9Na1qYacAAH8txIW5AkgO4zevToWq9u21BnnnlmnHnmmQUbDwAAgOIp6kWIAAAAYB8BFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkkpZhAYCI/OqZqocHAK2PHVAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQ6NPUEAKA2E/udeND7561dmmQeUAz1re8U/AwBTcEOKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhDqgUIv6aqM1h/ptANAYhXgtU0sUOFR2QAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJNQBBQCgQeqrJapOKLA/O6AAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEl0aOoJQGr1Fc0GAAojxWvuvLVLi/41gMKxAwoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJKEOKABQby1FNZRprvJZm2qFQvNhBxQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIIkOTT0BKDTF0gEKb97apQe933MvzVl967O+9Q0Ujh1QAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJdUABgEbLp46iWqE0V+qEQjp2QAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJNQBBaBZUncPaC7yqWHrOQvyYwcUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACCJDk09ASi0+gpB51NMGoohnyLlbWV9KtgOAG2THVAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAl1QAGaEfUxaanaSg1bABrHDigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQ6oABAvdT5hIOr72dEnWfYyw4oAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEh2aegKQ2ry1S+vtM7HfiUWfB61PPmsLmivPe9A4XgPSakvPWa1tbdkBBQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIokNTTwAAAGg7JvY7samn0KI09v9r3tqlBZlHodgBBQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSUAc0T8tf7qpmUT2aW42hxqjvsVgLbU9rWt8AUEx+T2pe8vl+pPw9J9kO6OLFi+Oaa66J8vLyKCsri9LS0jjiiCPik5/8ZEyePDl+/vOfx86dOw9pzPnz58fFF18cgwcPjq5du0ZZWVmccMIJccMNN8SKFSuK9EgAAABoiKLvgL7//vtx+eWXx8MPP3zAfVu3bo2tW7fGn//85/jFL34Rt9xyS/zsZz+L00477aBjfvjhhzFt2rSYO3dujdt37NgR77//fixbtix+8IMfxIwZM+LGG28s6OMBAACgYYoaQCsqKmLChAnx8ssvZ7f16tUrhg8fHsccc0xs2rQpXn/99XjzzTcjImLVqlUxYcKEWLBgQZxyyim1jllVVRWTJk2K+fPnZ7cNGzYsRowYERUVFbFw4cJYv359VFVVxU033RRVVVVxyy23FPNhAgAAkIeiBtC77rorC5/t2rWLb3/723H99dfHYYcdlvXJ5XIxd+7cuOqqq2LLli2xY8eOmDZtWvz3f/93rWPedtttWfjs3LlzzJw5My644ILs/l27dsU3v/nNuOuuuyIi4lvf+laMGTMmxowZU6yHCQAAQB6K+hnQmTNnZu1rr702vvGNb9QInxERJSUlccEFF8T/+3//L7vt1VdfjWXLlh0w3saNG+N73/tedvz973+/RviMiOjYsWPceeedMWXKlOw2b8MFAABoekULoB9++GG8/fbb2fGXv/zlg/Y/77zzokuXLtnxn/70pwP6zJ49O7Zv3x4REYMHD45p06bVOd6dd94Z7drtfXgvvPBCvPLKK4cyfQAAAAqsaAF027ZtNY6PPPLIg/Zv3759HHHEEdnxnj17Dujz6KOPZu2pU6dGSUlJneMNGDAgxo8fnx0/8sgj9cwYAACAYipaAO3Vq1d07tw5O3799dcP2n/jxo2xcePG7PjTn/50jfsrKytj8eLF2fHYsWPrnUP1PgsWLKi3PwAAAMVTtIsQlZaWxuc///ls5/G2226LiRMn1nibbXXTp0/Pdj3Hjx8fgwcPrnH/ypUrs/tLSkpi+PDh9c6hvLw8ay9fvrxBj4P8Nbcit8VU3+NQgLnlaS1rExrK8xrUzWvE33guaJ1q+74uz22OiJ0F/1pFvQrud77znXjiiSdi27Zt8corr8QJJ5wQN998c4waNSorw/Lqq6/Gv/7rv8Zzzz0XERFDhw6tcfGifVauXJm1e/fuXWN3tS4DBgzI2ps3b45NmzZFr169GvRYtsWW+GOuYbuoJ5WMa9B5AAAAjdHQDLMtthR4JnsVNYAef/zx8dxzz8U555wT77zzTqxatSqmTp1aa98jjzwyLrroovjOd75T47Og+7z33ntZu0+fPnl9/aOPPrrG8ebNmxscQHfH7tgSmxt0LgAAQFNobhmmqGVYIvZ+lvNPf/pT/PCHP4yuXbvW2W/ixIlx0UUX1Ro+I2pe1Gj/Ui512b/f/hdGAgAAIJ2i7oBGRGzatCm+9rWvxf333x9VVVVx9NFHx6hRo6JHjx6xZcuWePHFF+Ptt9+OuXPnxty5c2PatGnxox/9KNq3b19jnMrKyqzdsWPHvL52p06dahxXVFQ0/gEBAADQIEUNoG+88UacccYZ8e6770anTp3i7rvvjiuuuKJGuMzlcvGLX/wipk2bFlu2bImf/OQn0b59+/jRj35UY6zqn/nctWtXXl9/586aH5rNd+cUAACAwitaAP3oo49i0qRJ8e6770ZExE9+8pO4+OKLD+hXUlISkydPjp49e2Z1O+++++6YOnVqnHzyyVm/bt26Ze18dzL371d9jEPVPtpHt+je4PMBAABS6x5lDTpvW2yJ3bG7wLMpYgD95S9/Ga+99lpE7L0YUW3hs7px48bFhAkT4oknnoiIiJkzZ9YIoD169MjaGzZsyGsO69evr3FcVtaw//yIiG7R3dVsAQCAFqWhGeaPuQVFuYBR0QLo7373u6w9duzYvM4ZN25cFkBfeumlGvcNGTIka2/cuDEqKyvrLcWyevXqrF1WVtbgK+BSOI2tHdVS6nDlM091tNJqKWsHgPS8RvyN308otqJdBXffW28jau5eHkz1flu21Kw7M2TIkGjXbu90c7lcLF26tN7xlixZkrWHDh2a1xwAAAAojqIF0OoX/Nm8Ob+t2+q1Po888sga93Xu3DlGjhyZHT/99NP1jvfMM89k7XHjvH0WAACgKRUtgA4YMCBrP/XUU3mds2DBgqx93HHHHXD/eeedl7VnzZp10LHWrFkT8+fPr/VcAAAA0itaAD3zzDOz9ooVK+K+++47aP8FCxZkn/+MiJg4ceIBfS655JLo2rVrRESsXLky7r333jrH+9rXvha7d++9atOpp54a5eXlhzR/AAAACqtoAfQLX/hCjQsHTZs2LX784x9noXCfXC4XDz74YEyaNCm7rX///nHBBRccMGbv3r3j+uuvz46vvfbaePDBB2v02bVrV0yfPj3mzJmT3XbHHXc0+vEAAADQOEW7Cm6HDh1i9uzZMW7cuNixY0dUVlbG1VdfHd/+9rfjtNNOi549e8aWLVti8eLF8fbbb2fnderUKe6///7o1KlTrePefPPNsWjRoliwYEFUVFTElClT4vbbb4/y8vKorKyMhQsXxrp167L+M2bMiDFjxhTrYQIAAJCnogXQiIhTTjklnnrqqfhf/+t/xZ/+9KeIiFi3bl388pe/rLX/scceG/fdd1+MGjWqzjFLS0vj4YcfjmnTpmW7n8uWLYtly5Yd0O/WW2+Nm266qUCPBgAAgMYoagCNiDj55JPj9ddfj1//+tfx6KOPxksvvRRr166Nbdu2RdeuXaNPnz4xYsSIOPfcc+NLX/pSlJaW1jtm9+7dY+7cuXHFFVfE7Nmz44UXXoh169ZFaWlp9O/fPyZOnBiXXXaZ0is0a42tOdaW6nSpzwZAY3gdgeaj6AE0Yu/bcSdNmlTjc56FcOaZZ9a42BEAAADNV9EuQgQAAADVCaAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkESSOqBQKBP7nVhvn7ZSbLqtPE6geajvOSef52doKm1lfebzu0Fb+b+g+bIDCgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkoQ4orU599a3UzwQoPPUHoen5GaMlsAMKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkESHpp4AANA2zFu79KD3T+x3YpJ5ANB07IACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEmoAwoANAvqhAK0fnZAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJLo0NQTAADIx7y1Sxs9xsR+JzZ6DAAazg4oAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEOqAAQJuRTy1RtUIBiscOKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEh2aegIAAM3JvLVLm3oKBTGx34lNPQWAA9gBBQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJARQAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJLo0NQTgEKbt3ZpU08BAACohR1QAAAAkhBAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkOjT1BOBQzFu7tKmnAAAANJAdUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCXVAoYWa2O/EJF9H7VUAAArFDigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBIdmnoCtC3z1i5t6ik0GxP7ndjUU8hLY+fpew4AwD52QAEAAEhCAAUAACAJARQAAIAkBFAAAACSEEABAABIQgAFAAAgCQEUAACAJNQBpWDUewSA5qO+1+WWUo8aaF3sgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASagDmqctsTmezD2UHY8c0SkWPda/CWdEc5ayttofcwtiS2yucVv3KIuTSsYlm8PB5PN/oYZs6zTq7Hdi8cs7a9zmuZPmwNqkOWvur+vQWHZAAQAASEIABQAAIAkBFAAAgCQEUAAAAJIQQAEAAEhCAAUAACAJARQAAIAkBFAAAACS6NDUEwAAIL15a5fW22divxOLPg+gbbEDCgAAQBICKAAAAEkIoAAAACQhgAIAAJBESS6XyzX1JJqbPn36xMaNGw/ap1vXkhh2fMdEM6KlWf5y12Rfa1tsid2xu8Zt7aN9dIvuyebQWENHbG/qKVAEr63YFdu213yJ8dxJc2Bt5i/l6xl7tYbXdVqH2tZi7969Y8OGDY0aVwCtxeGHHx7btm1r6mkAAAA0G926dYutW7c2agxvwQUAACAJARQAAIAkBFAAAACS6NDUE2iOPvaxj8W7775b47YuXbrEoEGDmmhGAAAA6bz55puxY8eOGrd97GMfa/S4LkIEAABAEt6CCwAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYDWY9euXXHffffFWWedFQMHDozOnTtH375947TTTovvfve78de//rWpp0gr8fbbb8dPf/rT+MpXvhKf/vSn46ijjorS0tIoKyuLE044Ia688sp45plnGjT2/Pnz4+KLL47BgwdH165dszFvuOGGWLFiRYEfCW3JddddFyUlJdm/j3/844d0vrVJIS1ZsiSmT58en/nMZ6Jv377RqVOn6NevX5SXl8ell14a9913X6xfvz6vsaxNCmXx4sVxzTXXRHl5eZSVlUVpaWkcccQR8clPfjImT54cP//5z2Pnzp2HNKb1SYuWo07Lly/PDR8+PBcRdf7r3bt37vHHH2/qqdKCLVmyJHfyyScfdJ1V/zd27NjcX/7yl7zG3rJlS27KlCkHHa+0tDT3ne98p8iPktboxRdfzLVr167Geho4cGBe51qbFNKGDRtyF110UV7Poddcc81Bx7I2KZTNmzfnJk2alNe6/MQnPpFbtGhRvWNan7QGJblcLtfYENsarVmzJk455ZRYu3ZtRESUlJTE6NGj47jjjouNGzfGk08+GRUVFRERUVpaGr/97W9j/PjxTTllWqgHHnggvvzlL9e4bfDgwTFs2LDo2bNnfPDBB/H888/HmjVrsvv79esXzz77bAwaNKjOcauqquLzn/98zJ8/P7tt2LBhMWLEiKioqIiFCxfW2AmYMWNG3HLLLQV8ZLRmVVVVMWLEiFi2bFmN2wcOHBhvv/12vedamxTK6tWrY+zYsfHWW29ltx177LFRXl4ePXr0iIqKinjjjTdi6dKlUVlZGddcc038x3/8R61jWZsUSkVFRZx++unx8ssvZ7f16tUrhg8fHsccc0xs2rQpXn/99XjzzTez+7t06RILFiyIU045pdYxrU9ajaZOwM3V6NGja/xF/7//+79r3L9p06bc+PHjsz5lZWW5999/v2kmS4s2Z86cXETkjjvuuNy//uu/5tasWXNAn927d+fuvffeXJcuXbI1N3LkyNyePXvqHPfmm2/O+nbu3Dk3Z86cGvfv3Lkzd8MNN9T4q+nTTz9d8MdH63Tbbbdl6+bCCy88pB1Qa5NC+eCDD3KDBg3K1kl5eXmdu0hbt27NPfDAAwest+qsTQplxowZ2Rpp165d7vbbb8/t2LGjRp89e/bk5syZk+vevXvW94QTTqhzTOuT1kIArcXjjz+e/eB27Ngx9+qrr9bab9u2bTVe+G688cbEM6U1ePrpp3MzZ87MffTRR/X2ffjhh2u8sPzud7+rtd+GDRtyXbt2zfr9+Mc/rnPM6m/lOfXUUxv8OGg7li9fnuvUqVMuInIXXXRRbubMmXkHUGuTQrr88suzNTJ69Ojc9u3bGzyWtUkhffzjH8/WyP/+3//7oH0feuihGq/ttf3eaX3SmgigtTjrrLOyH9wrrrjioH3/8z//s8YuaFVVVaJZ0lZV/7zoP//zP9fa584778z6DB48+KA7pX/5y19qfI5vyZIlxZo6rcCePXtyo0aNykVE7qijjspt2LDhkAKotUmhvPLKK9naOPzww3PvvPNOo8azNimULVu21AiUL7744kH7f/TRRzXe4fTQQw8d0Mf6pDVxFdz9bNu2rcZ76//xH//xoP2/9KUvxeGHHx4REZs3b46FCxcWdX4watSorF3XZ+0effTRrD116tQoKSmpc7wBAwbU+PzyI4880ug50nrdfffdsWjRooiIuOuuu6J3796HdL61SaH8+Mc/ztqXXnppHHPMMY0az9qkULZt21bj+Mgjjzxo//bt28cRRxyRHe/Zs+eAPtYnrYkAup/nn38+uxR2165d46STTjpo/06dOsXIkSOz4wULFhR1flD9RWf37t0H3F9ZWRmLFy/OjseOHVvvmNX7WMPUZc2aNTF9+vSIiDj99NPj0ksvPaTzrU0KZffu3TFnzpzs+KKLLmrUeNYmhdSrV6/o3Llzdvz6668ftP/GjRtj48aN2fGnP/3pGvdbn7Q2Auh+li9fnrU/9alPRYcOHeo9p7y8vNbzoRiqX3W0f//+B9y/cuXK7K+nJSUlMXz48HrHtIbJx9VXXx1bt26Njh07xj333HPQv8DXxtqkUF577bX48MMPI2LvH4uHDx8eO3fujHvuuSfGjBkTvXv3js6dO8cxxxwTZ599dvz0pz+NXbt21TmetUkhlZaWxuc///ns+LbbbosdO3bU2X/69OnZ+hs/fnwMHjy4xv3WJ62NALqflStXZu2BAwfmdc6AAQOytuK/FNM777xT4y+ZZ5555gF9qq/hfb+E1af6Gt68eXNs2rSpkTOltXnggQfisccei4iIr3/96zF06NBDHsPapFD++Mc/Zu0hQ4bEqlWr4jOf+UxcddVVsXDhwti0aVPs3Lkz3n333Xj88cdj2rRpcfzxx8eSJUtqHc/apNC+853vRLdu3SIi4pVXXokTTjghZs+eHX/+85+jsrIy3nnnnXj88cfj9NNPj5kzZ0ZExNChQ7N2ddYnrU3923ttzHvvvZe1+/Tpk9c5Rx99dNbevHlzwecE+1x33XXZ224HDBgQ55xzzgF9GruGI/au4169ejViprQm7733Xnz1q1+NiIhPfvKT8Y1vfKPB4+xjbdIY77zzTtZu165dfPazn43Vq1dHRMTxxx8fJ510UrRv3z5effXVLHS+9dZbMXr06Hj22WcP2EGyNim0448/Pp577rk455xz4p133olVq1bF1KlTa+175JFHxkUXXRTf+c53anwWdB/rk9bGDuh+qn9w/LDDDsvrnOr99v/gORTK7Nmz45e//GV2fMcdd0SnTp0O6NfYNbz/GHDddddln0+65557al13+bA2KZQPPvgga7/00kuxevXqOOyww2Lu3LmxfPny+NnPfhYzZ86Ml19+ORYsWBA9e/aMiIjt27fHlClToqqqqsZ41ibF8OlPfzr+9Kc/xQ9/+MPo2rVrnf0mTpwYF110Ua3hM8L6pPURQPdTWVmZtTt27JjXOdV/GauoqCj4nOCll16Kq666KjueMmVKXHjhhbX2bewajrCO+Zvf//73cd9990VExCWXXBJnnHFGg8eyNimU7du3H3Db7NmzY/LkyQfcfsYZZ8Svf/3raNdu7688b7zxRtx///01+libFMOmTZvi6quvjuuvvz62b98eRx99dHzxi1+MadOmxZQpU+LjH/94RETMnTs3TjvttLjyyivrvLjgPtYnrYEAup/q76s/2AULqtt31dyI/P8yBfl666234pxzzslegD71qU/FPffcU2f/xq7hCOuYvbZv3x5XXnllRET06NEjvvvd7zZqPGuTQtn/M3AnnXRS/MM//EOd/U899dSYNGlSdvzAAw/UOZ61SSG88cYbMXz48Jg1a1a0a9cu7r777lizZk089NBDcc8998QDDzwQb775ZsydOze6d+8eERE/+clP4p//+Z8PGMv6pLURQPez7wPjEfn/tah6v+rnQ2OtW7cuJkyYEOvXr4+IiEGDBsW8efOyF6vaNHYN7z8Gbdc3vvGNrNbsv//7v2dvY2woa5NC2X8dnH/++fWeU73P888/X+d41iaN9dFHH8WkSZPi3XffjYi9wfKqq66K9u3b1+hXUlISkydPjocffji77e67744//OEPNfpZn7Q2Auh+evTokbU3bNiQ1zn7wkFERFlZWcHnRNv03nvvxYQJE2LVqlUREdG3b9948skno2/fvgc9r7FrOMI6JmLJkiXxwx/+MCL2voXxkksuafSY1iaFUn0tRUT83d/9Xb3nVO+zdevW2Lp1a63jWZs01i9/+ct47bXXImLvxYguvvjig/YfN25cTJgwITve/0q41ietjavg7mfIkCFZ+y9/+Ute5+y78l7E3icaaKwPP/wwPve5z2XFq3v06BFPPPFEHHvssfWeW30Nb9y4MSorK+u9ZHv1NVxWVuZKecSrr76a1Z1bvXp1jBw5ss6+1S/vv27duhp9b7755vjCF74QEdYmhbP/a20+uzv799m6dWscfvjhEWFtUli/+93vsvbYsWPzOmfcuHHxxBNPRMTe6z5UZ33S2gig+6le227ZsmXx0UcfRYcOB/9vql5XrCG18aC67du3x1lnnZW9AB1xxBExb968+Pu///u8zh8yZEi0a9cu9uzZE7lcLpYuXXrQ8BBhDXNwq1atynbi67Nr16548cUXs+Pq4dTapFCGDRtW47j6bmZd9u9T/aMM1iaFtO+ttxEH7tbXpXq/LVu21LjP+qS18Rbc/Zx22mnZlcO2b99+wF+h9rdz585YvHhxdjxu3Liizo/WrbKyMs4999xYtGhRRER06dIlfvOb38SIESPyHqNz5841Xpiefvrpes955plnsrY1TLFYmxTKscceG4MGDcqO/+d//qfec6r3KSsrq1EWw9qkkKpf8Cff+vDVa30eeeSRNe6zPmltBND9dOvWLcaPH58dz5o166D9H3744eyvqkcddVSMHj26mNOjFauqqoovfvGLsWDBgojYewn1X/3qVzFq1KhDHuu8887L2vWt4TVr1sT8+fNrPZe2a+rUqZHL5fL6V/3zSgMHDqxx3/6F161NCqX6RYUeeeSRevtX71Pba7W1SaEMGDAgaz/11FN5nbPvtT8i4rjjjjvgfuuTViXHAR577LFcROQiItexY8fca6+9Vmu/7du354477ris7/Tp0xPPlNbio48+yn3pS1/K1lKHDh1yv/rVrxo83oYNG3Jdu3bNxvvpT39aZ98vf/nLWb9TTz21wV+TtmvmzJnZGho4cOBB+1qbFMqf//znXGlpabZGHnzwwTr7Pv/887l27dplfR999NED+libFMqjjz6arY+IyP3sZz87aP/58+fX6D9r1qwD+liftCYCaB1OP/307If34x//eO7VV1+tcf9f//rX3IQJE7I+ZWVluffff79pJkuLtmfPntzUqVOztdSuXbvcnDlzGj3uzTffnI152GGH5ebOnVvj/p07d+a+/vWv13jRe/rppxv9dWl7DiWA5nLWJoXz1a9+NVsjXbp0qTWELliwINezZ8+s38iRI3N79uypdTxrk0KoqqrKDRkyJFsjnTt3zt199925jz76qEa/PXv25ObOnZvr3r171rd///65ysrKWse1PmktSnK5XK4xO6it1Zo1a+Lkk0+OdevWRUREu3btYsyYMTFo0KDYtGlTPPnkk7Fjx46IiOjQoUP87ne/q/HWXcjXj370o7jmmmuy409+8pPx2c9+Nq9ze/ToETNmzKj1vqqqqvjc5z5X4209n/rUp6K8vDwqKytj4cKF2fqOiJgxY0bccsstDXwUtGWzZs2Kf/zHf4yIvW/B3Vc7tC7WJoWyc+fOmDBhQjz77LPZbUOHDo2TTjop2rdvH6+++mq8/PLL2X19+/aNF198Mfr371/reNYmhfLiiy/GuHHjst8VI/auv9NOOy169uwZW7ZsicWLF9d4vuzUqVM88cQTcfrpp9c6pvVJq9HUCbg5W758ee7EE0+s8Zek/f/16tUr99hjjzX1VGnBvvWtbx10jR3sX327TR988EFu8uTJBx2jtLQ09y//8i9pHiyt0qHugOZy1iaF88EHH9R4y2Fd/0455ZTc6tWr8xrP2qQQXnzxxdzgwYPzej0/9thjc88991y9Y1qftAZ2QOuxa9eueOCBB2LOnDnx+uuvx4YNG+LII4+MQYMGxfnnnx+XXnpp9OzZs6mnSQt266231rmLWZ98dpsiIp588smYPXt2vPDCC7Fu3booLS2N/v37x8SJE+Oyyy5ziXYa5VB3QKuzNimUhQsXxs9+9rN47rnn4t13343du3dHnz59YuTIkTF58uQ477zzoqSkJO/xrE0K4aOPPopf//rX8eijj8ZLL70Ua9eujW3btkXXrl2jT58+MWLEiDj33HPjS1/6UpSWluY9rvVJSyaAAgAAkIQyLAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJCEAAoAAEASAigAAABJCKAAAAAkIYACAACQhAAKAABAEgIoAAAASQigAAAAJCGAAgAAkIQACgAAQBICKAAAAEkIoAAAACQhgAIAAJDE/wcFiJ8hSe58OAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 464,
       "width": 464
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(10)\n",
    "im = ps.generators.overlapping_spheres(shape=[100, 100, 100], r=10, porosity=0.5, maxiter=0)\n",
    "plt.imshow(im[:, :, 50]);\n",
    "im_poro = ps.metrics.porosity(im)\n",
    "print(f\"Porosity from image: {im_poro*100:.1f}%\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "snow = ps.networks.snow2(im, boundary_width = 0)\n",
    "network = snow.network\n",
    "pn = op.io.network_from_porespy(network)\n",
    "pn['pore.diameter']=network['pore.inscribed_diameter']\n",
    "pn['throat.diameter']=network['throat.inscribed_diameter']\n",
    "model=op.models.geometry.throat_length.cubes_and_cuboids\n",
    "pn.add_model(propname='throat.length',\n",
    "             model=model,\n",
    "             regen_mode='normal')\n",
    "model=op.models.geometry.pore_volume.cube\n",
    "pn.add_model(propname='pore.volume',\n",
    "             model=model,\n",
    "             regen_mode='normal')\n",
    "model=op.models.geometry.throat_volume.cuboid\n",
    "pn.add_model(propname='throat.volume',\n",
    "             model=model,\n",
    "             regen_mode='normal')\n",
    "pn.regenerate_models()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity from pnm: 67.5%\n"
     ]
    }
   ],
   "source": [
    "Vol_void = np.sum(pn['pore.volume'])+np.sum(pn['throat.volume'])\n",
    "Vol_bulk = 100**3 # from the image\n",
    "pnm_poro = Vol_void / Vol_bulk\n",
    "print(f\"Porosity from pnm: {pnm_poro*100:.1f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. Notes:: <Volume errors>\n",
    "\n",
    "\n",
    "    1) Note that in the example above, we assumed the original image was available for calculating the bulk volume. Assuming the original image is not available, calculating the bulk volume for the extracted network is another source of approximation error. Existing methods in topotools such as `get_domain_length` and `get_domain_area` tackle this problem to provide a better approximation of bulk volume of the network. \n",
    "\n",
    "    2) Note that in the example above, we assumed cubic pores and cuboid throats for the network. Assigning other diameters, shapes, and volume functions would result in a different estimation of void volume."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
